Model fluid in a porous medium: results for a Bethe lattice 
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We consider a lattice gas with quenched impurities or 'quenched-annealed binary mixture' on the 
Bethe lattice. The quenched part represents a porous matrix in which the (annealed) lattice gas 
resides. This model features the 3 main factors of fluids in random porous media: wetting, ran- 
domness and confinement. The recursive character of the Bethe lattice enables an exact treatment, 
whose key ingredient is an integral equation yielding the one-particle effective field distribution. Our 
analysis shows that this distribution consists of two essentially different parts. The first one is a 
continuous spectrum and corresponds to the macroscopic volume accessible to the fluid, the second 
is discrete and comes from finite closed cavities in the porous medium. Those closed cavities are in 
equilibrium with the bulk fluid within the grand canonical ensemble we use, but are inaccessible in 
real experimental situations. Fortunately, we are able to isolate their contributions. Separation of 
the discrete spectrum facilitates also the numerical solution of the main equation. The numerical 
calculations show that the continuous spectrum becomes more and more rough as the temperature 
decreases, and this limits the accuracy of the solution at low temperatures. 

PACS numbers: 05.50.+q, 64.70.Fx, 75.10.Nr 



I. INTRODUCTION 

When fluids are adsorbed in porous materials they be- 
have very differently from what we know in the bulk. 
This happens in both high-porosity materials, such as 
silica aerogels, and low-porosity materials, such as Vy- 
cor glass. Even an aerogel that occupies a few percent 
of the total volume significantly deforms and reduces the 
gas-liquid binodal of the fluid. At least three factors af- 
fect phase equilibrium of fluids in porous media: wetting, 
randomness, and confinement by the matrix. A number 
of models exist that emphasize one or two of those fac- 
tors, and they greatly facilitate understanding of possible 
behaviors of the fluid in the porous media (see review in 
Ref. [|). But theoretical models that take into account 
all three factors appear to be extremely hard to deal 
with. All existing microscopic theories are inconclusive 
even concerning the qualitative behavior of the gas-liquid 
binodal in the porous medium. Conventional approxima- 
tion schemes yield very different results depending on the 
model and on the level of approximation 0, [j| . Monte 
Carlo simulations suffer from long relaxation times and 
are performed in a very limited simulation box |lj, |J, |5( . 
This makes any exactly treatable model especially desir- 
able. We shall consider such a model, which features all 
the three factors. It derives from the well known lat- 
tice gas model of a fluid. We shall consider here only 
the simplest representation of a random porous medium 
(quenched impurities of equal size). The model is ex- 
actly solvable on the Bethe lattice, in the sense that 
we can derive a closed integral equation for the local 
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field distribution. Related approaches, leading to simi- 
lar equations, were previously known in the context of 
spin glasses 0, Q and the Random Field Ising model 
d, . The resulting integral equation has an interesting 
structure, which means that although numerical solution 
is relatively straightforward at high temperature, as tem- 
perature is lowered the field distribution (and hence also 
the numerics) becomes more and more involved. This 
problem is also worse close to the percolation threshold 
of the porous medium. 

The Bethe lattice is defined on a large uniformly 
branching tree of which the Bethe lattice is the part far 
from any perimeter site. (This is distinct from a Cayley 
tree which includes all parts.) The Bethe lattice should 
not be confused with the Bethe approximation, which 
usually denotes a set of mean-field-like self-consistency 
equations for order parameters. The Bethe approxima- 
tion is closely related to the cluster variation method or 
cluster approximation which derive these equations from 
truncated cluster expansion of either the free energy or 
entropy functional employing several variational param- 
eters. For many models these approximations become 
exact on the Bethe lattice, but not in the case studied 
here (except for special choices of parameters). This is a 
part of the motivation for the exact analysis presented be- 
low. Another important goal is to rectify a shortcoming 
of the grand canonical ensemble when dealing with fluids 
in porous media. This ensemble allows fluid to equili- 
brate throughout the pore space, including closed pores 
that are inaccessible to fluid in reality. For the Bethe 
lattice we are able to isolate and remove the effects of 
these closed pores. 

We shall formulate the model and present the basic 
equations in Section [H] The analysis of Section ITTT1 will 
bring forward the notion of finite clusters and their im- 
portance for a correct description of the fluid. Details 
of the numerical algorithm and computer-aided results 
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will be presented in Section llVl Section Ivl contains our 
conclusions. 



II. THE MODEL 

In lattice models of a fluid, the particles are allowed 
to occupy only those spatial positions which belong to 
sites of a chosen lattice. The configurational integral of a 
simple fluid is thereby replaced by the partition function 

Z = Trexp(-/?ft) , = l/(k B T), 

H = H - fj,N = --^2 !ij n i n j - ji^nu (1) 

ij i 

where rii, which equals or 1, is the number of particles 
at site i (i = 1 • • • V, where V is the number of the lattice 
sites [l0|). Tr means a summation over all occupation 
patterns. The total number N of particles is allowed 
to fluctuate; u is a chemical potential, which should be 
determined from the relation 

N=(j2 ni \ ;(■■•)„ = Z- 1 Tr(---)cM-m- (2) 

\ i I H 

Since particles cannot approach closer than the lattice 
spacing allows, lattice models automatically preserve 
one essential feature of the molecular interaction: non- 
overlapping of particles. The lattice fluid with nearest- 
neighbor attraction is known to demonstrate the gas- 
liquid transition only. Nevertheless, the lattice gas with 
interacting further neighbors possesses a realistic (which 
means argon-like) phase diagram with all the transi- 
tions between the gaseous, liquid, and solid phases being 
present . In this paper we deal only with fluid phases 
and restrict ourselves to the nearest neighbor interaction. 

The fluid adsorbed in the porous solid can be thought 
of as residing on a lattice, of which a fraction of sites are 
excluded or pre-occupied by particles of another sort. Al- 
though in practice these blocked sites, representing the 
solid, must form a connected network (if the solid is to 
remain static) we ignore this here and allow sites to be 
blocked at random. Thus, we have to consider the lat- 
tice gas with quenched impurities, whose Hamiltonian is 
given by 

U = -I ^ n i n j - K £ niX i "''E" 1 ' ( 3 ) 
{ij) {ij) 1 

where Yliij) means summation over all nearest neighbor 
sites, and quenched variables Xi describe the presence 
(x{ = 1) or absence (xi — 0) of a solid particle at site 
i. Each site of the lattice can be either empty (xi = 0, 
rii = 0), occupied by a fluid particle (x, = 0, n.; = 1), or 
filled with a particle belonging to the porous solid and 
therefore inaccessible to the fluid (xi = 1, rii = 0). The 
first term of Eq. J3J describes the nearest neighbor at- 
traction between the fluid particles (I > 0), the second 



one corresponds to the fluid-solid attraction (K > 0) or 
repulsion (K < 0) on the nearest neighbor sites. The 
hard-core repulsion is taken into account by mutual ex- 
clusiveness of particles: Xi + rii < 1- 

It is quite customary to establish the connection be- 
tween lattice gases and the Ising model. Indeed, a simple 
substitution Si — 2n; — 1 transforms into the Ising 
Hamiltonian. In the Ising model language Si = — 1 means 
an empty site, and Si — 1 corresponds to a site occupied 
by a fluid particle. The current model © is equivalent 
to the diluted Ising model with random surface field A: 

H = - ^ JijSiSj - ^2 hiSi + const, (4) 
{ij) 1 

Jij = JyiVf, hi = yiH + A;; A, = yjHy^Xj, (5) 

j 

J = 1 74; H = fi/2 + zI/A; H = K/2 - 7/4. (6) 

Here y,; = I — Xi describes the accessibility of site i to 
the fluid particles; the constant term in (JIJ embraces 
several pieces that depend only on quenched variables 
and therefore do not contribute to the thermodynamics; 
z is the coordination number of the lattice; the sum on 
j in (J5J spans nearest neighbor sites to site i. When 
deriving Q we used the equality rii = niyi. Aj is a 
random field correlated to surface sites: Aj ^ at the 
sites that are accessible to the fluid (yt = 1), and are in 
contact with the solid (at least one of the neighbor sites 
has Xj ^ 0). When H — 0, Eq. Q is a Hamiltonian 
of the diluted Ising ferromagnet, and we shall frequently 
refer to this limiting case throughout the paper. 

The grand thermodynamic potential of the model is 

F = -k B T{lnJ2---J2 e M-m)\ , (7) 

\ Si S v I {3.} 

1 1 

((•••)>{*}=£ ■•• E p({*i})(-)> ( 8 ) 

x-i— xv— 

where the distribution function p{{xi}) is a quenched one 
that describes the distribution of solid particles in space. 

We shall study this model on the Bethe lattice, which 
permits us to calculate an exact thermodynamic poten- 
tial. It is known that many problems on tree-like struc- 
tures, such as that depicted in Fig. ^ can be solved it- 
eratively. The term 'Bethe lattice' refers to an infinitely 
deep part of such a tree. 

Each site on the nth level has k = z — 1 neighbors at 
the (n+l)th level and one neighbor at the (n— l)th level. 
The partition function of the model can be calculated re- 
cursively, from upper levels to the bottom. For example, 
let us consider the cluster of sites framed by the dashed 
box in Fig.^ where k = 3. The following relation (which 
holds for all k) shows how tracing out spins at the upper 
level forms an 'effective field' on the lower level: 

k k 

£ ■ ■ ■ £ exp(/?S J2 JoiSi W = 

Si S k 1=1 1=0 
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FIG. 1: A schematic picture of the Bethe lattice with coordi- 
nation number z = 4. The arrowheads do not mean that the 
links are asymmetric, they just illustrate the way the partition 
function calculation proceeds (see text). 



[HC(J ol ,ht)}e X p(l3hoSo), 



1=1 



h = h Q + ^ U(Jqi, h, 



(9) 



(10) 



i=i 



U(J, h) = /3 _1 atanh[tanh(/3J) tanh(/3/i)], (11) 
cosh(/3J) cosh(/3h) 



C(J, h) = 2- 



cosh[/3J7(J, /i)] 



(12) 



We have thus obtained an effective field ho which is a 
sum of the original field ho and k contributions from the 
upper branches. We can proceed this way recursively 
downward, giving 



(13) 



at each level. In the non-random case, when hi = h and 
Jij = J, the effective field deep inside the tree satisfies 
the stable point relation that follows from Eq. (fH^ : 

h = h + kU(J,h). (14) 

For our model Q , we can expect that a properly defined 
effective field distribution tends to a stable limiting form 
at deep levels of the tree. 

Let us consider the probability Qi(y, h)dh that yi = y 
and hi = h: 

Q i (y,h) = (5 yyi S(h-k)) [ , (15) 

\ / {x} 

where 6(h) stands for the Dirac distribution, and 6 ab is a 
Kronecker's delta (6 a b equals 1 if a = b, and equals zero 
otherwise). Qi(0,h) is the effective field distribution at 
a matrix site, and Qi(l, h) is the field distribution for a 
site accessible to fluid. It is easy to see from (THfl) and (JjJJ 
that when yi = 0, Jij = and hi = 0, and the effective 
field hi equals zero, thus Qi(0,h) does not depend on i 
and has a simple form 



Qi(0,h) = Q(0,h) = co8(h) 
where we define 



(16) 



(17) 



so that Co is a fraction of the sites occupied by the ma- 
trix. When we move recursively from the surface of the 
tree deeper into its interior, Qi(l, h) changes, but should 
tend to a fixed point at infinitely deep levels of the tree: 
Qi(l, h) — ► Q(l,h). It follows from JT2J) that the result- 
ing distribution obeys 



Q(l, h) = Li [ [J] dh,Q( yj ,h 3 )]6(h -hi-J2 U(J t3 ,h 3 )) 
\ J i=i j=i 



(18) 



{*} 



In writing this equation, and throughout the following, 
we now specialize to the case where there is no correlation 
of the {yi} variables (which specify the sites occupied by 
the solid matrix) beyond nearest neighbor correlations, 
which are permitted. In this case, taking into account 
that k upper branches do not interact until they meet 
the chosen site, the fields {hj} are uncorrelated. (In more 
general cases, the yj at different branches are correlated; 
the joint distribution Q(yi, hi, • ■ • , yu, hu) cannot be de- 
coupled as a product of Q(yj,hj); and (|18fl is invalid.) 



Using the integral representation of the delta function 

6(h)=J^-eMm (19) 
we can factorize the delta function in l|18fl and get 

Q(l,h) - ( J ^ yi eM^(h-hi)) x (20) 

Y[ J dh J Q(y j , hj) exp(-i£U(Jij, h j))} 



j=i' 
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Configurational averaging with respect to the quenched 
variables (• • ■) s x \ in 1|2(J|) can now be performed explicitly, 
and it yields 

Q(h) = J ^Lexp(it(h-H))Q k (Z), (21) 

Q(0 = p exp(-itH)+p l J dh'Q(h')exp(-i^u(h')), 

where we introduced the notations Q(h) — Q(l,h)/ci; 
u(h) = U(J, h); p a — w a i/cr, and 

w ab = (5 ayi 6 by] ) {x} (22) 

is the probability to find yi = a and yj — b at a ran- 
domly chosen pair of nearest neighbor sites i and j. For 
example, Wu = (yiyj)i x \ is a probability that both sites 
in the pair are accessible to fluid. For each site accessi- 
ble to the fluid, p\ specifies the probability that a given 
nearest neighbor site is accessible as well; po = 1 — p\ is, 
correspondingly, the probability that this neighbor site is 
occupied by the solid. 

Note that Q(h) specifies the distribution of the effec- 
tive field created by k 'upper branches'. The complete 
one-site field deep inside the tree consists of contributions 
coming from all z nearest neighbors, and its distribution 
takes into account z equivalent branches 

Qz{l,h) = (yiSih-hi-^UiJij^hj))) 



= ci J |Uxp(^(/i-#))Q*(0, (23) 
Q z (0,h) = c Q 6{h). (24) 

This distribution determines the one-site average values. 
For example, in the magnetic interpretation of the model 
©, this total one-site field permits to calculate the av- 
erage magnetization (per magnetic site) 

m = ((ntSi^Jc! = J dhQ z (h)tanh(l3h). (25) 

Calculation of other thermodynamic properties is less 
straitforward. The thermodynamic potential of the tree 
is a sum of 'site' and 'link' contributions |6j,|7j: 

0F({x}) = fc^lnTr i exp(-/3H i )-^lnTr ij exp(-/W y -), 

< («) 

(26) 

where TLi and TLij are Hamiltonians of one site and a 
pair of sites, respectively, and the traces act only on the 
corresponding local degrees of freedom. The configura- 
tionally averaged thermodynamic potential of the Bethe 
latice must take into account only deep levels of the tree. 
We do this by assigning the stable point field distribu- 
tion to all the sites when doing configurational averaging. 
The Bethe lattice thermodynamic potential F is, there- 
fore, given by 



(3F/N = (3(F({x})) {x} /N 




where TLi and TLij depend on both fluid (S) and matrix 
(y or, equivalently, x) variables: 

TLi = hiSi', TLij = hiSi + hjSj + JijSiSj. (28) 

III. THE SOLUTION 

The result of the previous section is that thermody- 
namics of the model is determined by the one-site field 
distribution Q(h), and this distribution satisfies integral 
equation l|21(l . Equations of this type have been known 
previously in the context of spin glasses [|| Q and the 
Random Field fsing model [3, 01 ■ We did not encounter 
this kind of equation for the current model in the liter- 
ature, even for the limiting case of diluted ferromagnet 
(H = 0). In this section we investigate the ways of solv- 
ing this equation. 



(27) 

^ — — / dh i Q{y i ,h i ) I dhjQ(yj, hj) In Try exp(-/3W y ), 

y . =0 C Vi C Vj J J 

I 

For the diluted Ising ferromagnet in zero external field, 
when H = H = 0, Q(h) = 8(h) is always a solution of 
Eq. (|21|l . and it yields zero magnetization ((Si) = 0) 
or, in the fluid model language, the occupancy number 
equals one half ((rii) — 1/2) at all sites. This is known as 
the trivial solution in all calculations leading to an ana- 
lytical result, and this is the only solution at the temper- 
atures above critical. What are the nontrivial solutions? 

When there is no solid medium (c\ = 1) or when the 
solid sites are arranged in a nonporous block or slab with 
a smooth surface (pi = 1), the solution is Q(h) — 8(h — 
h), and h is given by Eq. (|14fl . This equation has a 
nontrivial solution (h ^ 0) at low temperatures: T < T c , 
T c = J/[/cBatanh(l/fc)]. In this case the results of the 
cluster (or Bethe) approximation become exact. 

When quenched chaos is present (pi ^ 1), and H = 
H — 0, the trivial solution becomes unstable at temper- 
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atures below 

T c = J/[fc B atanh(l/(fcpi))] (29) 

(this is an accurate result), signifying the appearance 
of the spontaneous magnetic order ((Si) ^ 0). Impor- 
tantly, though, the cluster approximation is not accurate 
when both magnetic ordering and quenched disorder are 
present. In the general case (H je 0, H / 0) we do 
not have any simple formula for the critical temperature 
or the field distribution. In what follows, we study the 
structure of this distribution more closely. 

First we shall consider the simple limitin g ca se of di- 
luted magnet without surface field (H = 0) [l^]- In this 
case matrix sites break the links between spins, but do 
not introduce any field that breaks the spin flip symme- 
try. It follows from (|29[1 that when p% < 1/fc, T c does 
not exist, and the spontaneous magnetic order does not 
appear. (The reason is that in this case the system is not 
percolated: there is no infinite cluster of linked spins, 
only finite clusters, and in finite systems the spin flip 
symmetry cannot be spontaneously broken.) Let us real- 
ize that in this system for any matrix (at any value of p\ 1 
except or 1) there is a finite fraction ^| of spins whose 
links are all broken, because all their neighboring sites 
are occupied by the matrix. The effective field at those 
sites equals just an external field H . Thus it is likely that 
the solution Q(h) contains a term in 8(h — H). In the 
case of zero external field (H = 0) one easily finds that 
a solution of the form Q(h) = Aq5(K) + q(h) does satisfy 
Eq. (|21|) . and one gets 

A = ( P o+PiA ) k , (30) 

— exp(i£h) x 

x[{p + PiA + q(Z)} k ~M> (31) 
= Pi J dtiq(ti)exp(-itu(h')), 

where we expect q(h) to be a nonsingular function, be- 
cause a delta function positioned in any other place 
(A x S(h — a), a ^ 0) or a set of such delta functions does 
not lead to selfconsistent equations of the required form. 

Let us show that Aq has a simple physical mean- 
ing, namely, it equals the probability that all fc upper 
branches are finite. Consider first the probability Pf that 
a branch is finite. It consists of two possibilities: either 
the first link is broken (with probability po), or it is un- 
broken, but all farther branches connected to it are finite 
(with probability piPj). This yields Pf = po + piP k . 
Then, the probability that fc given branches are finite 
equals A^ = P^ and satisfies Eq. (|30|1 . Obviously, 
Eq. (|30|l always has the solution Ao — 1, which leads 
to the trivial solution for Q(h) because of the normaliza- 
tion condition J Q(h)dh = 1. Ao = 1 is the only solution 
when pi < 1/fc, as discussed above. At the percolation 
point (p\ = 1/fc) it is a two-fold solution, which signifies 
a bifurcation. When pi > 1/fc, a solution for Aq in the 
interval (0,1) always exists. 



When H je 0, the previous ansatz that Q(h) — A8(h — 
H) + q(h) does not work: 5(h — H) placed into the rhs of 
Eq. i|21|) generates, under iteration of the equation, a set 
of k different delta functions, as does any delta function 
positioned in some other place. This seemingly rules out 
the idea of finding the discrete levels in this case, but a 
fully solvable case z = 2 (the linear chain: see Appendix 
|A"|) resurrects the hope. In that case solution contains an 
infinite series of delta functions. Therefore, in the general 
case we shall seek the solution in a form of a sum of an 
infinite set of delta functions and a nonsingular function: 

oo 

Q(h) = Y,ai5(h-hi) + q(h). (32) 
i=i 

We now show that such a solution can really be obtained. 
The point is that, similarly to the zero field case (H = 0), 
the equation for the discrete spectrum separates from the 
equation for the non-singular part q(h). Placing l|32l) into 
Eq. J2U, we get the equation for the discrete spectrum 

jr ai 6(h-hi) = f^exp^(h-H)) x (33) 
i=i J 71 

oo 

[po +pi^ai cxp(-i£,u(hi))] k . 
i=i 

The number of delta functions is infinite, but their cu- 
mulative weight remains finite, and 

oo 

J2 a l= A °- ( 34 ) 
i=i 

The latter relation can be obtained by integrating 
Eq. (|33|l with respect to h and noting that the total 
weight of delta functions satisfies Eq. Ij30|l . Eq. Ij34(l 
means that the obtained discrete spectrum is just the sin- 
gle h = level that we had in the case H = split by the 
external field. Each (a;, hi) pair corresponds to a certain 
finite tree 'growing' up from the given site. For example, 
(a\,h\) — (pq,H) corresponds to the 'zero' tree, when all 
fc links to the upper sites are broken; naturally hi = H , 
and ai equals the probability that all the neighboring fc 
sites at the upper level are inaccessible. 

Now let us turn back to the most general case, when 
both the surface and the external fields are present (H ^ 
0, H ^ 0). Inserting into (|2TJ 

OO 

Q(h) = s(h) + q(h), s(h) = aiS(h - hi) (35) 

i=i 

one can obtain explicit expressions for the discrete s(h) 
and continuous q(h) parts of the spectrum: 

s(h) = J ^expmh-H))§ k (H), (36) 
q(h) = J |Uxp(i£(/i-ff)) x (37) 

x[{m + m} k -s k (a 
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where 



HO = Po exp(-i£ff) +px ai exp(-i£u(hi)), (38) 



z=i 



9(6 = Pi dh'q(h')exp(-i£u(h')). 



(39) 



Note that Eq. (|38(l for the discrete spectrum has a 
closed form and does not depend on q(h), whereas the 
equation for q(h) does include s(h). This is because the 
discrete spectrum comes from finite trees, whereas the 
continuous spectrum is an attribute of the infinite vol- 
ume. Finite branches know nothing about the infinite 
volume, therefore s(h) does not depend on q(h). At the 
same time, at any given site finite branches may connect 
onto the infinite tree; therefore the equation for q{h) does 
depend on the discrete spectrum s{h). 

Let us consider the structure of Eq. 1)3 7|l. It can be 
rewritten in the form 



q(h) 



k 

f gexp(i£(/i-ff))£ 
J n 1=1 



ZiO^iO- (40) 



Let us recall that by construction of equation (|2~TJl for 
Q(h) if we insert some field distribution as Q(h) in the 
rhs , we obtain in the Ihs the field distribution one level 
deeper the tree. The same property holds for Eq. (|40|l . 
Then Eq. (|40|l permits a clear interpretation: each term 
of mixes into the next-level distribution the possibil- 
ity that / branches are infinite, with the other k — I finite; 
I runs from 1 to k, the I — term was rightly eliminated 
by the last term in Eq. 1)37(1 , because this term would give 
rise to the possibility of a finite tree and would give a dis- 
crete spectrum contribution. This interpretation makes 
s(£) and <?(£) responsible for finite and infinite branches 
respectively. The surface field H does not enter Eq. I|37|) 
explicitly; its influence on q(h) is mediated by the dis- 
crete spectrum s{h). In a way, s(h) (or, equivalently, 
s(0 UM) encapsulates the surface influence. The only 
explicit reminder of a porous medium in the equation for 
the continuous spectrum is the p\ multiplier in Eq. (|39|1 . 
indicating that links between sites are not always present. 

The complete effective field acting on an accessible site 
Q z (l,h) (defined by ((23(1 ) likewise consists of a discrete 
spectrum s z (h), and a continuous one q z {h): 



Q z (l,h) = c x s z {h) + c x q z {h) 



2tt 



exp(i£(/i — H)) x 



(41) 
(42) 

(43) 



xM0 + mr-s z (a 

The total weight of the discrete part 



s z {h)dh = A z /k (44) 

equals the probability that all z branches connected to 
this site are finite. When this happens, the site belongs 
to a finite closed volume. 

When considering the fluid in the grand canonical en- 
semble (as we have done so far) all sites are in equi- 
librium with an implicit bulk fluid, whose temperature 
and chemical potential are parameters that describe this 
equilibrium. But in real experiments fluid particles can- 
not access any disconnected pore volumes (described by 
s z (h)), and those parts of the system do not respond to 
changes of the chemical potential. Were we considering a 
magnetic system, spins in those finite separated clusters 
would be able to react to changes of external field H, and 
thermodynamic potential 127(1 would be valid. In con- 
trast, in order to describe the experimental situation for 
a fluid in porous medium, we have to exclude the contri- 
bution of finite volumes to the thermodynamic potential. 
Fortunately, we are now able to separate this contribu- 
tion because finite branches always generate a discrete 
spectrum, in contrast to the infinite cluster, which yields 
the continuous distribution. 

Inserting 1(35(1 and l(41|) into the expression for the ther- 
modynamic potential 1(27(1 , we can identify three distinct 
contributions of finite volumes 



feci / dhiS z (hi) In Tr 4 exp(-^7ij) 



y<=i> 3/3=1 



zwio / dhis{hi) lnTr^ exp(— j3H.i 



ion / dhis(hi) / dhjsQij) lriTty exp(— jSTLij) 



The first term above corresponds to an accessible site 
whose branches are all finite, the second one describes a 
linked pair of accessible sites with outer branches finite, 
and in the third term one site of the pair is accessible to 
the fluid, but all its branches are finite. Note that the 
division into these three terms does not reflect the pres- 
ence of different types of site but reflects the subdivision 
into site and link contributions in 1(261) . 

Dropping these 'finite volume' contributions (and con- 
stant terms corresponding to the sites occupied by the 
matrix) we obtain the thermodynamic potential of the 
fluid in the accessible volume as: 



0F/N = feci / dfnqMlnTnexpi-pHi 
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"2 Wn I dhiq(hi) j dhj[q(hj) + 2s(hj)] lnTr^- exp(-f3Hij) 



-zwio I dhiq(hi)]xLTrijexp(—f3Hij) 

I 



(45) 



The procedure of 'contribution classification' is espe- 
cially simple when we consider the expression for the den- 
sity (per free volume) n = ((nA^) i x \/c\. The grand 
canonical ensemble (GCE) predicts 



1 fri _j_ 

dh[q z {h) + Sz (h)]--[l + t^M(3h)\ = — Z 



1 



(46) 



whereas when we take into account that closed finite 
pores are inaccessible to the fluid the microscopic fluid 
density at those sites has to be set to zero, and we obtain 
the corrected, 'infinite cluster' (IC) expression 



= J dh{q z (h) • i[l + tanh(/3ft)] + s z {h) • 0} 

= [1-1 + J dhq z (h)tanh(/3h)}/2 
m + 1 



(47) 



An, 



where x is a fraction of free sites belonging to the finite 
cavities, defined by Eq. I|44|) . and the difference between 
the GCE and IC densities 



An= dhs z (h) • -[1 + tanh(/3ft)] 



(48) 



is simply the number of fluid particles in finite cavities 
per number of the free sites. 



IV. NUMERICAL CALCULATIONS 

In order to solve the main equation 121fl we have to 
calculate the discrete spectrum (|36H and use it when 
solving for the continuous part of the field distribition 
via l)37|) . The procedure of the discrete spectrum calcu- 
lation is iterative. We place the zeroth approximation 
s(ft) = p%5(h - H - kH) into the rhs of Eq. $® and 
obtain in the Ihs a sum of k + 1 delta functions, one 
of which coincides with the initial one and has the cor- 
rect weight Pq. Placing the new approximation into the 
rhs of Eq. ll.'Uill recursively, we can obtain as many terms 
as we need. The unpleasant feature is that the number 
of levels in the spectrum grows exponentially during the 
iterations, and many levels have negligible weight. For- 
tunately, one can drop the low-weight levels at each step, 
because they produce only lower-weight levels during the 
following iterations. 

Let us note that the weights of delta functions depend 
only on parameter p\ , which is determined by the struc- 
ture of the solid matrix. The positions of the delta func- 
tions depend also on other parameters (T, H, H, and 
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FIG. 2: The coordinates and weights of delta peaks in the 
dicrete spectrum s(h) for k = 2, H = 0, pi = 0.9, T/T c = 0.8; 
Aq = 1/81; the 132 heaviest levels are shown for H/J — 0.1, 
H/J = 0.4, and H/J — 1. The weights for all the 3 cases are 
the same, only the positions of the delta peaks differ. 



J) . As we discussed in the previous section, the discrete 
spectrum corresponds to finite trees growing from a site. 
Each level in the spectrum corresponds to a tree of a cer- 
tain size and form, and its weight is a probability to find 
such a tree in the system. Fig.[21shows an example of the 
discrete spectrum for the case of a relatively sparse solid 
that blocks only a small fraction of sites. In that case 
the 20 'heaviest' levels carry more than 99.96% of the 
total weight. As p\ approaches the percolation thresh- 
old, the convergence of the series becomes progressively 
poorer. For example, the 200 heaviest levels in the case 
k = 2, pi = 0.6, contain only 87% of the discrete spec- 
trum's weight. This naturally agrees with the fact that 
there appear a lot of large finite trees near the percolation 
point, and the number of forms of a tree quickly (com- 
binatorially) grows with its size. One other interesting 
feature of the discrete spectrum is its bandlike structure, 
with a\ densely grouped around certain values, and hi 
subdividing into several branches as H increases. 

Once s(h) is calculated to a given high accuracy (de- 
fined by e = (Aq — J^i °z)A4o), we can proceed by solving 
equation (|37|l . A numerical algorithm must represent the 
function q(h) as an M- vector q(hi), hi = h 1TL in + «A/i, 
Aft = (ft max - ft m in)/M, i = • • -M - 1, with M suffi- 
ciently large. It is helpful to note that q(h) is non-zero 
only on a finite interval [ft m i n , ft max ]: this permits use of 
Fourier series instead of Fourier integrals. The lower and 
upper limits of this interval may be estimated as 



min(-fcj, (k- 1)H — J), 




FIG. 3: The continuous part q(h) of the effective field distri- 
bution for a weakly dilute Ising ferromagnet (co = 0.1) with 
H = 0, at different temperatures below the critical point. 



max(fc<7, (k - 1)H + J). 



(49) 



These bounds are far from exact at high temperatures, 
but this does not make a significant impact on the effi- 
ciency of the calculations [l5j . 

The numerical solution of Eq. 1)3 7|l is an iterative pro- 
cess. We start from some initial distribution (this may 
be a uniform distribution) and insert it into the rhs of 
Eq. JsSZJ. The resulting Ihs is a new approximation to 
the solution. Each iteration produces the field distri- 
bution q(h) one level deeper the tree. Therefore such 
iterations must converge provided the stable point dis- 
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FIG. 4: The continuous part q(h) of the effective field distri- 
bution for a strongly dilute (but still percolated) Ising ferro- 
magnet (co = 0.4) in zero external held at different tempera- 
tures below the critical point. Some pictures contain a dashed 
line, which shows the h > 1 part of q(h) at an enlarged scale 
(plotted against the right y axes). 



tribution defined by Eq. 121|l exists and numerical errors 
introduced by the discretization are small enough. Naive 
quadratures for Eq. (|21|l would result in a formula 



M-l 

l(hi) = exp^/i, -H)) x 

j = l-M 

M-l 
1=0 

27Tj/(/i max ^min); 



(50) 
(51) 



which takes of order M 3 calculations at each iteration, 
because each of the indices i, j, and I in (|50() and l|51|) 
take M essentially distinct values. This quadrature is 
both inefficient and incorrect. First, to implement the 
correct numerical algorythm, it is important to note that 
exp(— i^u(h')) is a quickly oscillating function of h! when 
£ is large, and this requires a special quadrature when cal- 
culating the h! integral in l|37|) . Second, the Fast Fourier 
Transform algorithm takes only 0(M In M) steps to per- 
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FIG. 5: The high-density branch of the GCE binodal (n+) for 
the model with H = at different matrix densities cq. The 
low density branch is symmetric: n_ = 1 — n+. The perco- 
lation concentration for the lattice used (k — 2) is Co = 1/2. 
The inset shows the dependence of n + at zero temperature 
on the matrix density. The dotted lines represent the corre- 
sponding results of the cluster variation method [l7j|. 



form the summation over j in l|5U|) for all i. Also, Press 
et al 0| describe an algorithm that permits simultane- 
ous calculation of M integrals over b! in 1|37|) using again 
only 0(M In M) operations. Therefore, with advanced 
algorithms, each iteration costs 0(M In M) arithmetic 
operations. This is a dramatic gain with respect to the 
quadratures estimate of 0(M 3 ), because as seen below, 
M is typically many thousand. 

The numerical results we report here are limited to 
coordination number 3 (k = 2), and a completely random 
matrix, when there is no correlation between quenched 
variables even at the nearest neighbor sites: hence w a t = 
c a Cb and p a — c a - Since we are studying the behavior of 
a fluid in a porous media, we can set aside very dense 
matrices with p\ < 1/fc; these have no infinite connected 
volume accessible to the fluid, so that q(h) = 0, and 
thermodynamic potential (|45|l equals zero. 

First we shall consider the simpler case of a diluted 
Ising ferromagnet (H = 0), which corresponds to fluid in 
the porous matrix with a moderate fluid-matrix attrac- 
tion (K = 1/2). In this case the exact critical tempera- 
ture is known from Eq. J2HJ); when H = and T > T c , 
Q(h) — 5(h) and m = (n = 1/2), while at lower 
temperatures there are two symmetric non-trivial solu- 
tions Q+(h) = Q-(-h), with Q+(h) = A 5(h) + q+(h), 
Ao = (co/ci) 2 . The temperature dependence of q+(h) 
is shown in Fig. EI One can see that the q+(h) line re- 
peats itself, but not in a trivial fashion. The continuous 
field distribution acquires even more structure for denser 
matrices (Fig.0J. At low temperatures q+(h) becomes a 
set of increasingly sharp peaks, and this greatly increases 
the errors introduced by the discretization in the numer- 
ical algorithm. Therefore this algorithm does not work 
at very low temperatures. Nevertheless, using extremely 



fine discretization (we used M up to 2 17 = 131072), we 
covered an interesting range of temperatures. With the 
algorithm just described we were able to reach T/T c = 
0.4. In the magnetic interpretation of the model © the 
solutions Q+(h) and correspond to exactly op- 

posite magnetizations (m+ = — m_), which reflects the 
fact that direction of the spontaneous (zero-field) magne- 
tization is not predetermined. In the fluid interpretation 
these two solutions form the two branches of the binodal 
(GCE: ?i + = (to+ + 1)/2, n_ = l-n+). Fig.|31shows the 
resulting GCE binodals at different densities of the ma- 
trix. According to Eq. (|48|l . the 'infinite cluster' binodals 
for this particular case (H = 0) have the same size and 
shape as GCE binodals, they are only shifted to lower 
densities by An = at/2, An = 6.8 ■ 10~ 4 , 7.8 ■ 10" 3 , 0.039, 
0.15 for Co = 0.1, 0.2, 0.3, and 0.4, respectively. In the 
general case (H ^ 0), the density shift between the GCE 
and IC binodals depends on temperature. 

The binodal density (or zero field magnetization, in the 
Ising language) at zero temperature can be found from 
the fact that in this case all the spins in the infinite clus- 
ter are aligned, whereas finite clusters do not contribute 
to the spontaneous magnetization, and the magnetiza- 
tion (per magnetic site) equals plus or minus the frac- 
tion of magnetic sites belonging to the infinite cluster: 
m + (T = 0) = 1 — x. This gives the resulting half-width 
of the binodal, which is depicted in the inset of Fig.JSJas a 
function of matrix density cq. One can see that the clus- 
ter approximation produces the binodals that are flatter 
at the top and have an incorrect width. For a lattice with 
only sparsely distributed solid sites, the cluster approxi- 
mation results are quite accurate. 

Now, let us consider instead the situation where there 
is no matrix- fluid attraction (hard core matrix): K = 0, 
H = —J. Fig. shows how the field distribution at 
constant chemical potential (fixed H) now changes with 
temperature. In the case depicted there is a competition 
between the surface field H = — J < and the bulk field 
H = 0.21J > 0. At high temperatures the effective field 
prefers the sign of the bulk field and is mainly positive. 
At T/T c < 0.82, an alternative solution for q(h), which 
favors negative fields, corresponds to the global minimum 
of the thermodynamic potential and becomes preferable. 
At T/T c = 0.82 the two solutions yield the same value 
of thermodynamic potential at different fluid densities, 
and correspond to the two branches of the binodal. It 
should be noted that T c is defined as given by Eq. I|29() 
and is not at (but is close to) the critical temperature for 
this matrix. Fig. shows the temperature dependence 
of the fluid density at different values of H. The jumps 
correspond to coexistent densities at the binodal. One 
can see that in this case the binodal is not symmetric. 
The results depicted in Fig. [7| correspond to the 'grand 
canonical fluid' described by Eqs. H27fl . I|46|l . The fraction 
of the free sites that belong to finite closed cavities (x — 
(co/ci) 3 ) is as small as 0.14% for this relatively sparse 
matrix, therefore the binodal of the 'infinite cluster fluid', 
which excludes their contribution, is indistinguishable on 
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FIG. 6: The effective field distribution Q(h) for a fluid in a 
sparse hard-core matrix (H = — J, Co = 0.1) at H = 0.21J 
at different temperatures. The crosses show weights ai and 
locations hi of leading delta functions in the discrete part of 
spectrum s(h) = aiS(h — ft;), plotted against the left y 
axes. The lines depict the continuous part of the distribution 
q(h) (against the right y axes). The x axis is h — H. 




FIG. 7: The temperature dependence of the equilibrium fluid 
density at different values of H for the model with H — — J 
at matrix density Co = 0.1. The numbers at the dashed lines 
are the values of the external field H . The lines jump over 
the phase separation region, delimited by the solid line. The 
binodal within the cluster approximation fl^l is depicted with 
a dotted line. 
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FIG. 8: The field distribution for the model with H = —J 
at matrix density c = 0.3, H = 0.61, and T/T c = 0.64 
calculated with truncated discrete spectrum s(h). One of the 
cases depicted (a) corresponds to 48 leading delta-functions 
taken into account in s(h) (all the terms with weights a; > 
10 -4 ). In the other case (b) there are 804 terms in s(h) (the 
terms with weights a ; > 10 -6 ) 



the plot from that depicted by the whole line in Fig. 
Let us also note that due to a special symmetry of the 
model 01 i the binodal of the symmetrical case H = J 
(K = 2 J, attractive matrix) coinsides with that in Fig.[7| 
flipped with respect to n = 1/2 line. 

Next, we would like to consider denser matrices, and 
here the poor convergence of s(h) leads to a problem. Wc 
did not have this problem in the case of dilute magnet 
(H = 0), because its entire binodal, in view of spin- flip 
symmetry of that special case, corresponded to H = 0, 
when the discrete spectrum consisted of only one term: 
s(h) — Ao5(h), which could be exactly calculated for 
the matrix of any density (Eq. J5DJ). For the hard core 
matrix (H — — J) this no longer takes place. For Co = 0.3 
the 48 leading delta-functions of s{h) contain 95.5% of its 
weight, 804 terms contain 98.6%. If we simply ignore the 
residual weight and solve Eq. (|37fl with an incomplete 
s(h), the resulting q(h) strongly depends on the number 
of terms taken into account in s(h) (see Fig. |SJ). 

To overcome the problem, we use the following trick. 
First we calculate the leading 21 terms of s(h). They 
correspond to 21 relatively small finite trees, because the 
larger trees have smaller weight. Then we determine the I 
lightest terms among them and add A = (Aq — Y^iLi a i)/l 
to their weights. The resulting discrete distribution 



s{h) = a.iS(h - hi) 



i=l 



21 

J2(ai+&)S(h-hi), (52) 

i=l+l 



has the correct total weight. This way we distribute the 
missing weight of large trees among the largest of those 
taken into account. From the physical point of view this 
means that we replace the 'very large' branches by the 
'large' ones and in this way forbid 'very wide' pores. The 
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FIG. 9: The same as Fig. |H| but calculated with an adjusted 
discrete spectrum s(/i), Eq. I."i2t . The arrows point to the 
groups of ai's raised by A. The resulting continuous spectra 
corresponding to (a) 48 and (b) 804 terms in s(h) are indis- 
tinguishable. 



latter implies a certain correlation in the matrix and is 
always an approximation, since when deriving equation 
(|37|l we permitted only the nearest neighbor matrix cor- 
relations. Anyway, the results must become accurate in 
the limit I — ► oo, and we find that the method converges 
already at small I (Fig. |5J). A remarkable new feature 
of q(h) in Fig. is the presence of steps, accompanied 
by sharp peaks, especially visible at h — H = —1.85 or 
— 1.75. It is possible that these peaks are not accurate 
(see Appendix[BJ) , they could be numerical artefacts aris- 
ing in places where q(h) in reality has only a step |20| . 
Note anyway that these peaks are not delta functions: 
their heights remain almost unchanged when we increase 
resolution from M = 2 14 to M — 2 15 , whereas if they 
were delta peaks, their heights would then double pi). 

Having overcome this convergence problem with s(h), 
we can calculate phase diagrams for denser matrices 
(Tigs. HUIllfl . For these denser matrices, unlike the case 
Co = 0.1, the difference between the grand canonical en- 
semble binodals and the infinite cluster ones becomes 
clearly visible, although the fraction of finite closed cav- 
ities is still not large (x(co — 0.2) = 1.6%, x(cq = 0.3) = 
7.9%). This fraction grows to 100% at the percolation 
point (which happens at cq = 0.5 for this lattice), thus 
making a big difference near percolation. 

One can see that GCE and IC binodals for H = —J 
are indiscernable, which indicates that for this case closed 
pores in GCE are almost empty, and their contribution is 
negligible [2^|. On the contrary, for H = J, they are al- 
most saturated with fluid, and An is close to the fraction 
x of closed pores in the matrix. This situation should 
also hold for strongly attractive or repulsive matrices: 
\H\ > 1. For dense enough matrices with smaller \H\, 
An will be nonzero for both negative and positive H and 
will depend on temperature. 

In the above results for our numerical treatment of 
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FIG. 10: The binodals for the models with H — — J and 
H = J at matrix density Co = 0.2. The bold lines (c) repre- 
sent the 'infinite cluster' binodals. The thin dashed lines (b) 
are the grand canonical ensemble binodals (in the H = — J 
case this line is hidden by the overlapping (c) line). The 
binodals within the cluster approximation flq| are depicted 
with a dotted line (a). The pairs of (a) and (b) binodals are 
symmetric with respect to n = 1/2. The pair of (c) lines is 
symmetric with respect to n = (1 — x)/2 = 0.492. The small 
squares on (b) and (c) lines are the actual data points through 
which the interpolating lines were drawn. 
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FIG. 11: The same as Fig. 1101 at matrix density co = 0.3. 
The pair of (c) binodals is symmetric with respect to n = 
(1 - x)/2 = 0.461. 



the exact integral equations, deviations from the clus- 
ter approximation are clearly visible in the phase dia- 
grams (Figs. 0| , and become very large in some cases 
fFigs. HUI lllfl . Note that the cluster approximation per- 
mits no simple way round the failure of the grand canon- 
ical ensemble to prevent fluid from entering closed pores. 

It is interesting to note that some other related ap- 
proaches 0, predict an unusual double-hump binodal 
of the fluid in moderately dense strongly attractive or 
repulsive matrices (in terms of the current model, this 
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means large Kierlick et al |3J considered an off- 

lattice molecular model of fluid in a quenched disordered 
configuration of spheres on the basis of a replica symmet- 
ric Ornstein-Zernike equation. They studied a sequence 
of approximations, the first of which corresponded to the 
mean spherical approximation (MSA), and predicted a 
binodal of the usual form. The higher approximation 
yielded a binodal with two critical points or even two 
disconnected binodals. Their highest order approxima- 
tion again gave only one critical point, but predicted a 
shoulder on one of the binodal's sides. The MSA applied 
to the model J2J) on the five-dimensional fee or bee lattice 
again predicted a double-humped binodal. A problem 
identified for MSA was that for the 3-dimensional lattices 
it predicted a double-humped binodal even for the bulk 
fluid without any porous matrix. Therefore the authors 
°f could not conclude definitely about the existence 
of two critical points for these models, despite the simi- 
larity of their results to the double-humped binodals seen 
in Monte-Carlo simulations on similar models 0, Q . Our 
theory does not give any double- humped binodal, at least 
in the range of parameters we have studied. 

V. CONCLUSIONS 

We considered the Bethe lattice model of fluid in a 
porous medium. The recursive character of the Bethe 
lattice permits an exact treatment, whose key ingredient 
is an integral equation for the effective field distribution. 
Solutions of this equation consist of a sum of discrete 
and continuous spectra, and these spectra have distinct 
physical interpretations. The discrete spectrum comes 
from disconnected finite pore spaces, whereas the con- 
tinuous spectrum is a contribution of the infinite pore 
space which in reality is the only one accessible to fluid. 
The continuous spectrum develops more and more struc- 
ture at low temperature, which means that a numerical 
solution for it becomes impractical below a certain tem- 
perature. However, the physical results found for tem- 



where v(h) is the function inverse to u(h), that is 
v(u(hj) = h, and v(h) is its derivative. The above 
equation says that Q(h) equals poS(h — H — H) for h 
in the close vicinity of H + H, otherwise one has to 
look at the value of Q at v(h — H) and multiply it by 
P\v{h — H). This observation leads to a solution of the 
form Q(h) — J^i <^i5{h — hi). The leading delta function 
is given explicitly in the rhs of Eq. (|A1|) . and positions 



peratures above this threshold are both consistent and 
reasonable. Despite use of a Bethe lattice they differ 
significantly from results calculated using the cluster (or 
Bethe) approximation, which cannot handle the complex- 
ity of the field distributions that we find. 

The dichotomy between the two types of pore space 
mentioned above is not exclusive to the Bethe lattice, but 
universal. Any microscopic model of fluid in a random 
porous media that uses the grand canonical ensemble will 
include contributions of the finite cavities, unless this is 
carefully subtracted off, as we managed to do for the 
Bethe lattice. In the grand canonical ensemble, these 
cavities are in equilibrium with the external bulk fluid, 
but in real-world experiments they are inaccessible and 
do not respond to changes of chemical potential. This 
marks an important distinction between models of fluids 
in porous media and disordered magnetic model to which 
they are equivalent in the grand canonical ensemble; for 
magnetism, finite clusters do contribute to the free energy 
and it is right to include them. 
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APPENDIX A: SOLUTION OF THE INTEGRAL 
EQUATION FOR THE FIELD DISTRIBUTION 
FOR THE ONE-DIMENSIONAL CHAIN 

In the one-dimensional chain k = 1, and Eq. (|21|) can 
be solved analytically. Performing first the integration 
with respect to £ in the rhs of Eq. (|21|l , one gets 



(Al) 



of the further delta functions are given by relations 

v{h 2 -H) = h i; v{h 3 -H) = h 2 ; ••• (A2) 

This means, in particular, that in the vicinity of h 2 

Q(h) = piv{h - H)Q(v(h - H)) 

= Pl v(h - H) Po S(v(h -H)- ht) 

= pip S(h - h 2 ), (A3) 



Q(h) = Po S(h -H-H)+ Pl J dh'Q(h')S(h -H- u(h')) 

= PQ 6{h — H — H) + Pl v(h - H)Q(v(h - if)), 

I 
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FIG. 12: A function with discontinuities (periodic on 
[—0.5; 0.5]) and its Fourier expansion, truncated at M — 16 
and M = 32. 



and therefore = PiPo- In general, one can find that 
ai = pia/_i. Finally, we can write the solution in the 
form 



Q(h)= Po J2p[~ 1 5(h-hi) 



(A4) 



i=i 



where hi are given by the recursion that follows from 
Eq. E3, 

hi=u[hi-i)+H;l>2;hi = H + H. (A5) 



APPENDIX B: ON THE NATURE OF SPIKES IN 
THE CONTINUOUS DISTRIBUTION q(h) 



In our numerical calculation we represent q(h) as an 
M- vector and use truncated Fourier series. This is guar- 
antied to work well only for smooth functions. For exam- 
ple, truncated Fourier series result in peaks near disconti- 
nuities (Fig. El- These are called 'Gibbs ears'. Note that 
Gibbs ears have a height that remains fixed, while their 
width narrows to zero as the resolution improves. There 
is accordingly no area beneath a Gibbs ear, unlike a delta 
function. The spikes in Fig. ED do not look exactly like 
Gibbs ears: Gibbs ears appear in pairs, and the positive 
ear has a negative counterpart, whereas spikes in Fig. El 
are asymmetric, with negative spikes much smaller and 
almost absent at the fine discretization (M = 2 16 ) that 
we used to produce the plot. It is still possible that equa- 
tion/recursion H37f) somehow creates a type of positive- 
only Gibbs ear effect, but this does not seem likely 
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